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Abstract. A stable approach for integrating the impedance matrix in cylindrical, radial 
inhomogeneous structures is developed and studied. A Stroh-like system using the time- 
harmonic displacement-traction state vector is used to derive the Riccati matrix differential 
equation involving the impedance matrix. It is found that the resulting equation is stiff and 
leads to exponential instabilities. A stable scheme for integration is found in which a local 
expansion is performed by combining the matricant and impedance matrices. This method 
offers a stable solution for fully anisotropic materials, which was previously unavailable in the 
literature. Several approximation schemes for integrating the Riccati equation in cylindrical 
coordinates are considered: exponential, Magnus, Taylor series, Lagrange polynomials, with 
numerical examples indicating that the exponential scheme performs best. The impedance 
matrix is compared with solutions involving Buchwald potentials in which the material is 
limited to piecewise constant transverse isotropy. Lastly a scattering example is considered 
and compared with the literature. 



1. Introduction 

Wave propagation in layered elastic media has been widely studied resulting in a variety of 
solution approaches. These include the use of scalar and vector potentials [T], the transfer 
matrix method [21 EH HI El IE] , an d the delta matrix method [8] . Alternative, computationally 
stable methods have also been developed, e.g. the stiffness matrix [HI HO], the global matrix 
[TT] . and the reflectivity method [12]. Such approaches are limited to isotropic or transversely 
isotropic materials whereas we are interested in general anisotropic solids in order to develop 
scattering solutions related to metamaterial devices such as acoustic cloaks [131 El US] which 
can be modeled as radially inhomogeneous anisotropic solids. The goal of this paper is to 
produce a stable solution method for such materials. 

We considered a matricant propagator solution [161 [IT] based on Stroh formalism to solve 
for scattering from a generally anisotropic material. The method involves the creation of a 
state space representation of the system where six first order, ordinary differential equations, 
must be integrated which may often diverge. A stable scheme for finding the solution to 
the differential equations, inspired by [H], is developed by combining the matricant and 
impedance matrices. The scheme is considered with several different expansion methods to 
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yield relatively high orders of accuracy and is compared with solutions from the conditional 
impedance matrix and Buchwald's scalar potentials |19j . 

The outline of the paper is as follows. In section EJ we begin with definitions of the impedance 
and matricant matrices. An explicit method for finding the impedance in piecewise uniform, 
transversely isotropic materials is developed in section [3j This method also serves as a tool to 
compare with more general solution methods based on the Riccati matrix differential equation 
for the impedance matrix. It is found that the Riccati differential equation is stiff and leads 
to exponentially growing instabilities. In section HI an alternative approach for integrating 
the matricant is derived to find a stable means of integrating the impedance matrix. We then 
use different expansion techniques used in the integration process in order to find higher order 
accurate schemes. Lastly in section [5] a scattering example from the literature is considered. 



We consider time harmonic wave motion in radially inhomogeneous cylindrically anisotropic 
solids. The associated equilibrium equations for linear elastodynamics in cylindrical coordi- 
nates are summarized in|A] Here we need only focus on the relation between the vectors U(r) 
and V(r) associated with displacement and traction, respectively. Precise definitions follow 
from Eqs. ( Il6l) (which includes the superscript n that is here omitted for simplicity). The 
dimension of each vector is taken as m, where m is either 3, 2 or 1; m = 3 in general, m = 2 
if z— dependence is not considered, and m = 1 for pure out-of-plane shear horizontal (SH) 
motion. For the moment we may consider m as general. The main focus of the paper is the 
m x m conditional impedance matrix z defined such that 



It is shown in Eq. (117)1 that the equations for linear elastodynamics can be cast as a system 
of 2m linear ordinary differential equations [TT] 



where Q + = — TQT, + denotes the Hermitian transpose and T is defined in (15T1) . It follows 
from Eqs. ()T|) and (j5J) that z(r) satisfies a differential Riccati equation 

dz 

h zQi — Q 4 z — zzQ 2 z — zQ 3 = 0, (3) 

dr 

with assumed initial condition z(ro) at some specified r = ro, hence the name conditional 
impedance. 



2. Impedance and matricant matrices 






(2) 



STABLE IMPEDANCE METHODS 



3 



One approach to solving for the conditional impedance matrix, z, is to first solve for the 
2m x 2m matricant M(r, r ) which is defined as the solution of the initial value problem 

^(r,r ) = Q(r)M(r,r ), M(r ,r ) = I (2m) , M = ^ . (4) 

Hence 

rj(r) = M(r,r )?7(r ). (5) 

Using the relations 

U (r) = (Mi - iM 2 z(r )) U (r„) , V (r) = (M 3 - iM 4 z(r )) U (r ) , (6) 

which follow from ([TJ), the conditional impedance can be expressed in terms of the matricant 

as 

z(r) = i(Mg - M 4 z(r )) (M 2 - iM 2 z(r )) _1 . (7) 
The propagator nature of the matricant is apparent from Eq. (J5j) and from the property 
M(r, ri) M(ri,ro) = M(r, ro), and in particular M(r, ro) = M(ro, r)^ 1 . Also, the symmetry 
(l5Tj) implies M(r,r ) = TM + (r , r)T. Hence, M _1 (r,r ) = TM+(r, r )T, that is, M is 
T-unitary [20]. 

An alternative approach to finding z uses the two point impedance matrix, which by defi- 
nition relates the traction and displacement vectors at two values of r according to [17] 

-V(r)J -- lZ ^ r °^U(r)J' Z "^Z 3 Z 4 

The two point impedance matrix has the important property that it is Hermitian, Z = Z + 
[TTj . The relations between the matricant of (J4j) and the impedance matrix of ([8]) evaluated 
at cylindrical surfaces r, r are easily deduced [T7] 



M(r,r ; 



Z(r, r ; 



— Z 2 1 Z 1 1Z2 1 
iZ 3 — iZ 4 Z 2 Zi — Z 4 Z 2 y 

-iM^ 1 ]^! iM^ 1 



M4M2 M] - M 3 — iM 4 M2~ 
Introducing ([9]) into ([7]), we can relate the conditional impedance z(r) to the two point 
impedance matrix Z(r, r ) according to 



z(r) = Z 3 (Z 1 -z(r )) % - Z 4 . (10) 
3. PlECEWISE UNIFORM TRANSVERSE ISOTROPY 

In this section we develop an approach suitable for piecewise uniform transversely isotropic 
cylinders by explicit calculation of the global impedance matrix Z of (jHJ), from which the 
conditional impedance can be found using fTTOl) . The method is based on a recursive algorithm 
proposed by Rokhlin et al. [9] called the stiffness matrix method. The analysis in [9] was 
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restricted to multilayered media in Cartesian coordinates, whereas the present method is 
applicable to cylindrically layered media of transverse isotropy. We will refer to Rokhlin and 
Wang [9] several times in this section to note the similarities and differences of the approaches. 



Z 





Figure 1. A cylindrically anisotropic multilayered media is considered in the 
system of cylindrical coordinates. The media consists of n anisotropic layers 
with different densities and elasticity tensors in general. 

Consider n > 1 layers of uniform transversely isotropic materials with the kth layer rk-i < 



r < rk, k G l,n, see Figure [U The explicit form of the local two point impedance matrix of 
the fcth layer is Z h (rk,rk-i) as defined by Eq. ( l60l) . Denote the global two point impedance 
matrix for the cylinder between r and by Z K = Z (rfc,ro). Our objective is the global 
two point impedance matrix for the entire cylinder, Z(r n ,ro) = Z N (r n ,ro). 

Consider first the two bordering layers between r = tq and r = r 2 and sharing the r = n 
surface. Continuity of displacements and traction on the interface implies 



V \ (Zf Z§\ U 

-vj - [z* zi) lu a 



V x \ _ . Z\ Z\\ /Ux 



-v 2 i - *\Z b 3 Z\) vu 2 



(11a) 
(lib) 



where Z a = Z 1 (ri,ro), Z b = Z 2 (r2,ri). From the second row of Eq. filial) and the first row of 
Eq. flllbj) . we have 

Ux = - {Z\ + Z\Y l (Z§ U + Z\ U 2 ) . (12) 
Introducing Eq. ffl2|) into Eqs. filial) and flllbl) . we define the impedance matrix Z 2 (r2,ro) 
that relates the traction vector to the displacement vector on the inner (r = ro) and outer 
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r 2 ) surfaces of the bilayer, 



-vj=- iz >- r »)(u:i 



where 



(rja Tifrja i r 7b\ ^rja r 7a( r 7a i '7b\ ^rjb 

-z|(zs + z?) _1 zg z\ - z\(z% + z\y x T 

Z a i = Zj, Z\ = Z\ and Z\ (k = 1, 2, z = T~4) are given by Eqs. QBQ). Note that Eqs. ( lTTa|) and 
( lllbj) are similar, apart from a sign change, to Eqs. (19) and (20) in [9]. 

Employing (1T31) recursively, the global impedance matrix Z K (rfc, ro) for the cylinder between 
ro and rfc is obtained with 3x3 components 

(rvK-l yK-l/yit, yK-n-^K-l 7K-I frvfc , 7K-1\ - 1 rjk \ 

~ ^2 {^1 + ^4 J ^3 _Zj 2 V^l + ^4 J ^2 / 1 r\ 

-z^ + zf^zf 1 z\-z\{z\ + zf- i y l z k 2 J ' 1 J 

where Z^-" 1 , (i = 1, 4) are the 3x3 sub-matrices of the matrix Z K_1 (rfc_i, ro) for k — 1 layers, 
Z*, (i = 1,4) are the 3x3 sub-matrices of the matrix Z fc (rfc,rfc_i) for the fcth layer, defined 
by Eq. (|60|) . The global impedance matrix for the iV-layered cylinder is obtained by using Eq. 
(USD (N - 1) times. 

The main differences between the present results and those of [9] are, first that by con- 
struction the local Z k and global Z K two point impedance matrices are Hermitian matrices. 
Secondly, the present results are valid for cylindrically layered structures, as compared with 
those of [9] which are for multilayered structures in Cartesian coordinates. Despite the differ- 
ences, we note that the two point impedance matrix Z k of Eq. (jSJ) and the global two point 
impedance matrix Z K (rfc,r ) are, apart from some sign differences, similar to the stiffness 
matrix K m and global stiffness matrix K M of Rokhlin and Wang [9]. 

4. Stable solution technique for general anisotropy 

In this section we propose and demonstrate a stable numerical scheme for solving for the 
conditional impedance matrix z(r) defined by (CQ) in the case of arbitrary radially dependent 
cylindrical anisotropy, i.e. density and elastic moduli are arbitrary functions of r: p(r), Cijki{r). 

4.1. Stability issues. Direct numerical integration of either Q for the conditional impedance 
z(r) or (J4]) for the matricant M(r, r ) is not a feasible strategy. The stiff nature of (J4]) leads 
to exponentially growing instabilities for M. These become unavoidable at large values of n 
and/or kr for the elastic problem. On the other hand, singularities and numeric instabilities 
may form when Eq. ([3]) is integrated, which is a well known issue for Riccati equations [18] . 
The singularities (poles) of the impedance matrix occur at finite values of r associated with 
traction-free modes for the given frequency. One can, in principle, avoid the singularity by 
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switching to the differential equation for the inverse of the impedance: the admittance A = z 1 
[2T| p. 136]. The admittance satisfies a Riccati differential equation, which follows from Eq. 
(U as 



dA 

dr 



+ iAQ 3 A + AQ 4 - Qi A + iQ 2 = 0. 



(16) 



However singularities again arise, this time corresponding to rigid modes. One could develop 
a numerical scheme that switches back and forth between integrating the impedance and 
admittance but since the locations of the singularities are not known in advance, and in 
practice they can be very close together, this does not offer a viable method. The curves 
shown in Figure [2] exemplify the problem of singularities in the impedance. The appearance 
of the peaks in Figure [2] indicate values of r beyond which accurate numerical solution cannot 
be obtained regardless of the step size in the integration scheme. The inability of such standard 
methods to obtain correct results was the motivation behind the proposed solution technique 
discussed next. 
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Figure 2. Solid, aluminium cylinder integrated with 200 evenly spaced steps, 
using the fourth order scheme of section 14.4.2} from r = 0.5 to r = 1.0, with 
k z = and ka = 10. Plotted is the determinant of the upper left 2x2 sub- 
matrix of the the 3x3 conditional impedance matrix normalized by n 3 + 1 vs. 
r. Where Eq. (|58|) was used for the initial impedance at r = 0.5. 



4.2. Mobius scheme. Numerical difficulties in integrating the matrix Riccati equations are 
well known and have been studied extensively. A procedure mentioned earlier for avoiding 
singularities is a generalization of the idea of switching, where one undergoes a change of 
variables, which is closely related to invariant embedding methods, see e.g. [22]. Here we 
follow a different approach, based on [18], which views the solution of the Riccati equation as 
a " Grassmanian flow" of m- dimensional subspaces (the impedance matrix) on a larger vector 
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space of dimension 2m. The idea is to recast the equation for z in the form of a forward 
marching scheme of step size h based on 

z(r + h)= i(M 3 - iM 4 z(r)) (Mi - iM a z(r))~\ (17) 

where M = M(r + h, r). The key to the method is that M can always be calculated in a 
stable manner for sufficiently small step size h. This approach is one of a class of methods 
called Mobius schemes, which by design are formulated on the natural geometrical setting 
of the larger vector space, in this case that of M. Accordingly Mobius schemes are able to 
handle numerical instability and pass smoothly and accurately through singularities |18j. The 
method therefore combines both the matricant and the impedance, each of which is unstable 
when solved in a global sense separately. 

4.3. Approximations for M(r + h, r). The Mobius scheme shifts the problem to that of 
finding approximations for M(r + h, r) accurate to some given order in h. The Peano series 
[20] of cascading integrals for the matricant is formally guaranteed to remain bounded for any 
h, but is numerically impractical. Our objective is to develop approximations for M(r + h, r) 
in the form 

M(r + h,r) = I (2m) + hM^(r) + ^'(r) + ^M^(r) + . . . , (18) 

where the terms M^ 2 -*(r) do not require explicit integration schemes for their evaluation. 
Consider first the case of Q sufficiently smooth. Then the identity 

ft 

M(r + h,r) = I (2m) + J Q(r + s)M(r + s, r) d s, (19) 

o 

may be written in the form of a series in powers of h by using (|18p for the left member and a 
Taylor series evaluated at r for Q in the integral, 

ft 



hM W + |W 2 ) + |W 3 ) + . . . = I (JQ + S Q' + ^Q" 

o 

+ fjcr + ...] x [i (2m) + S M« + ^+...]jd S) 



(20) 



with the argument r understood for all functions. Comparing equal orders of h k in fTSOj) implies 
M (1) (r) = Q(r), 
M (2) (r) = Q'(r) + Q(r)M (1) (r), 
M (3) (r) = Q"(r) + 2Q'(r)M (1) (r) + Q(r)M (2) (r), 



M (4) (r) = Q'"(r) + 3Q"(r)M (1) (r) + 3Q'(r)M (2) (r) + Q(r)M (3) ( 



r 
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etc. The matricant may now be found by employing ffl8l) . The series ff2TT) . while illustrative, 
is restricted to profiles that are analytically smooth functions of r. It is not suitable for 
piece- wise constant or piece- wise smooth profiles, which are of practical interest. Derivatives 
of the profile are therefore to be avoided if possible. In that sense, the approximation formed 
from Eqs. ( |T8|) and (}2~Tj) is only valid to 0(h), and the iterative scheme (fTTI) shares the same 
accuracy. 

An expansion accurate to second order can be obtained by using a Taylor series expansion 
evaluated at the midpoint [18] 



Q(r + s) « Q(r + ^) + (s- ^)Q> + + + 0(s 3 ). (22) 

Substitution into (TT9]1 then gives, using f[P8"j) . 

M (1) = M< 2 > = Q(r + £), = (I (2m) + iq'(r + ~))M«. (23) 

This leads to an expansion up to 0(h 2 ) requiring only Q at a single position with no derivatives 

h s h 2 ^ 9 . h , 

The form of (|24p suggests an alternative expression that is accurate to the same order 



TS2nd: M(r + h, r) = I (2m) + hQ(r + -) + — Q 2 (r + -) + 0(h 3 ). (24) 



EXP 2nd(a) : M(r + h, r) = exp (hQ(r + |)) + 0(/i 3 ). (25) 

Approximations (|24|) and (|25|) . together with Eq. (fT7|) each yield a second order accurate 
Mobius scheme. EXP 2 nd has the added feature that it is unimodular, and hence energy 
conserving. Detailed comparisons are provided in section 14.61 



4.4. Lagrange interpolation expansions. We now consider Lagrange polynomial expan- 
sions for Q in Eq. (fl9|) in order to obtain higher order expressions without using derivatives 
of Q(r). The Lagrange polynomial of order n approximates Q(r + s) with the expression 



n n / \ 

Q(r + s) = J2Q(r + x J h)L J £), L s (x) = ]J l^-), (26) 

3=0 11 1=0,1^3 KXj %lJ 

where Xj G [0,1], j = 0,1,..., n are chosen points. Substituting into fl20l) and using the 
notation = Q(r + Xjh) implies the sequence 

M ( fc ) = | Lf ] Q Xj | M^" 1 ), MM = I, 

k = 1,2,3.... (27) 

Lf ] =k I L i (x)x fe - 1 dx, 
Jo 

Note that ^" =0 -^j = 1- I n the following subsections we derive expansions based on fl26l) for 
n = 1 and n = 3. 



STABLE IMPEDANCE METHODS 



9 



4.4.1. Two point approximation: Halves. We approximate Q with two points using ( 126]) for 
n = l. In this case the integrals Lj can be simplified with the result that 




Taking equi-space points xq — -g and X\ — |, yields 

h h 2 

LP 2nd: M(r + h, r) = I (2m) + - (Qi + Qs ) + — (Qi + 5Qa) (Qi + Qa ) + 0(h 3 ). (29) 

jr x 4 4: Za^X- ^ 4 4 4 

This again gives a Mobius scheme of second order accuracy in h. It also suggests, by analogy 
with ((25D, 

EXP 2nd(6) : M(r + h, r) = exp (^Qs) exp (^Qi) + 0(/i 3 ). (30) 

2 4 2 4 

Note that the expansions fl29l) and fl30|) only agree with one another to first order, 0(h). 



4.4.2. Four point approximation: (Fourths). Taking four evenly spaced points to approximate 
Q, Xj = -g+ 4> j — 0, 1, 2, 3, and using the symbolic algebra program Maple, yields 

M(1) =lQ^ + §Qf + §Qt + iQ^ 



M( 3 )=(-^Q | + iQ | + XQ | + ilQ | )M( 2 ), 

M (4) =( - 1Q| + HQ f - tIoQ! + iliQl)M (3) . 

Substitution of these terms into f|T8l) gives M(r + /i, r) up to fourth order accuracy. Inter- 
estingly, when more points were taken to evaluate Q the numerical accuracy was not found 
to improve. This was tried with even spacings, using from five up to ten points. Finally, by 
analogy with (130 j) we define 

EXP 2nd(c) : M(r + h,r) = exp (|Qr) exp (f Q|) exp (|Q|) exp (|Qi) + 0(h 3 ), (32) 

which, like ( l29l) and ( J30l) . is consistent with the four-term Lagrange interpolation scheme ( l3lT) 
to first order. 



4.5. Fourth order Magnus integrator scheme. The Magnus integrator, created by Wil- 
helm Magnus [28], and further developed in [26] with a convergence proof and recurrence 
relations, is a method to approximate a solution to (j4j) with 

M(r + h,r) = e n M(r,r-h). (33) 
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Here we consider a fourth order Magnus integrator scheme similarly done in [27] for the 
Helmholtz equation. We use the following definitions to march a solution forward in r, 

M(r + h,r) = e"M(r, r - h), 

MG4th: n = ^(Q(i) + Q(2)) + -^-(Q(2)Q(i)-Q(i)Q(2)), (34) 

Q(D = Q(r + h{\ - ^)), Q (2 ) = Q(r + h{\ + ^)). 

As a fourth order scheme the numerical precision of this method is very similar to that of 
the four point Lagrange polynomial approximation (I3ip . which is seen in the examples of the 
following section. 




Figure 3. Solid, aluminium cylinder integrated with 2000 evenly spaced steps 
from r = 0.5 to r = 1.0, with k z = 0, n = and ka = 5. Plotted is the difference 
of the determinants of the upper left 2x2 sub-matrices of Eq. fl58l) and those 
calculated from section HI As noted in section 14.61 the right hand side refers to 
the type and order accuracy, they are listed top to bottom as worst to best at 
r — 1. 

4.6. Numerical examples and convergence. In order to illustrate the convergence prop- 
erty of the different expansions proposed (exponential, Magnus, Taylor series, Lagrange poly- 
nomials), we consider a solid aluminium sample with properties p = 2.7kg /m 3 , E = 70GPa, 
G = 2QGPa and radius of r = 1 and normalized these properties with respect to water for 
which p = 1.0kg/m 3 and speed of sound in water c = lA70km/ s. The numerical values re- 
ported were computed by implementing the Mobius scheme 017|) and in the case of the Magnus 
integrator implementing ([7]), starting at r = 0.5 with initial condition given by the explicit 
solution fl58|) from [T7J and discussed in §|4j In all examples we take k z = 0. Figure [3] plots 
the difference of the determinant of the upper left 2x2 sub-matrix of the exact conditional 
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impedance from (15 8|) . and that calculated by iterating f lT7|) until r = 1 is reached. The right 
hand side of Figure |3] refers to the type of approximation used: Taylor Series (TS), Lagrange 
Polynomial (LP), Exponential (EXP), Magnus (MG), and to the accuracy order. Thus, LP 
1st was calculated using Eq. (EL}, TS 2nd by $23}, EXP (a) by $25}, EXP (b) by (ETiD, EXP 
(c) by ([32]), LP 2nd by ([29]), LP 4th by (13"!]) . MG 4th by ([34]), and LP 3rd was calculated 
using a Lagrange Polynomial with points {xo, X\, X2} — {g, |, §}■ Interestingly, the three EXP 
methods (Eqs. ( 1251) . ( l30i) and ( 1321) ) gave similar results and were the best results for the fewest 
number of approximation points. Figure [4] plots the difference of the upper 2x2 sub-matrices 
at r = 1.0 vs. the number of steps used in the iteration from r = 0.5 to r = 1.0. As expected, 
the higher order schemes are more accurate and require fewer steps in the integration process 
to yield the same accuracy as a lower order scheme. 




3500 4000 



# of steps 



Figure 4. Solid, aluminium cylinder integrated from r = 0.5 to r = 1.0, with 
k z = 0, n = and ka = 25. Plotted is the difference of the determinants of the 
upper left 2x2 sub-matrices of Eq. (1551) and those calculated from section @] 
at r = 1 vs. # of steps. 



5. Scattering Example 

In this section we explore the use of the impedance matrix by considering acoustic scattering 
from a solid aluminium cylinder immersed in water. Perpendicular plane wave incidence, 
i.e. k z = 0, in a uniform exterior fluid is considered with time harmonic dependence e _KJ * 
henceforth omitted. The total radial stress and displacement fields in the surrounding fluid 
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Figure 5. Plotted pressure field described by Eq. (|35|) for an aluminium cylin- 
der, integrated from r = 0.5 to r = 1.0 (area between the two red circles drawn) 
using the fourth order scheme, Eq. (}3~Tj) . with 500 steps. The initial impedance 
at r = 0.5 was found using Eq. (1581) . ka = 5, k z = 0, and a to t = 2.468. 



are 



a rr = -Kk J2 i n {Ukr) + B n H£\kr))e 

11= — OO 

OO 

«r = - E i n (X(^) + J B„^ 1) '(A;r))e i ^, 



(35) 



where r is the radial coordinate, K is the bulk modulus, k is the wave number, Hn^ is the 
Hankel function of the first kind, and the coefficients B n are to be determined pQ . We use the 
definition of the conditional impedance matrix, noted in Eq. ([!]), and write this statement for 
the innermost radial coordinate at r = b, for which we find the initial impedance matrix from 
Eq. (158]) . and for the outer surface at r = a where 

V(6) = -izxU(6), V(a) = -iz 2 U(a). (36) 

The conditional impedance matrix, z 2 = z(a), is found from the integration technique outlined 
in section HI or for transversely isotropic materials may be found directly using Eq. fl58|) . 
Considering acoustic fluid in the exterior we write Eq. (|36l) for r = a in detail for which the 
shear stress, a r g, must be zero with 



ui 



-iz 2 



Pi <?1 



j \P2 <?2 

Eliminating ug using the second row of Eq. (1371) implies 



<j rr i . ^ 



-IZq. 



(37) 



(3f 
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Using (155)) and equating with (138)) we find the scattering coefficient 

= KkaJ n (ka) + z j' n (ka) 
KkaH { n ] (ka) + z Hi iy (ka) ' 
Numerical simulation was conducted for a solid aluminium cylinder with properties normalized 
with respect to water and with the total pressure field illustrated in Figure [51 Figure [6] shows 
both the total scattering cross-section a tot and the back-scattering amplitude f(ir), where the 
far field form function, f(6) is 



f(9) = ^=Y, i2n ~ le " R " cosn6 > 

* n=0 



where eo = 1 and e m = 2 for m > 0. The total scattering cross section is then 

4:71 

°tot = t~ Imag(/(0)). 
ka 

Figure [6] closely matches the behavior of a similar analysis conducted in [23] . 



(40) 



(41) 



ka 

Figure 6. Total scattering cross section and backscattering amplitude plot- 
ted against non-dimensional frequency, ka. Aluminum cylinder with the same 
properties as listed in figure integrated using the fourth order scheme from 
r = 0.5 to r = 1.0 using 500 steps. The initial impedance at r = 0.5 was found 
using Eq. (j5SJ) . 



6. Conclusion 

Two computationally stable methods were considered to calculate the impedance matrix. 
Seeking solutions of 3D elasticity in the form of time-harmonic cylindrical waves, a matrix 
Riccati equation for the impedance matrix was formulated. Direct integration of the Riccati 
equation is numerically unstable, it leads to exponentially growing instabilities, which is in- 
escapable at high frequency, large values of n, and layer thickness. We integrated the Riccati 
equation for the impedance matrix which has singularities along the real radial coordinate 
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associated with traction free modes. To avoid instabilities we developed a new stable numer- 
ical scheme for cylindrically anisotropic structures that passes through these singularities by 
combining the impedance and matricant. This scheme evaluates the impedance matrix for 
continuous systems by integrating the Riccati equation over the thickness of each layer. Differ- 
ent expansion methods were considered and compared, it was noted that matrix exponential 
approximations yielded the best results. 

An alternative method was developed to obtain the global impedance matrix for anisotropic, 
cylindrically layered media using the impedance matrix for each layer. The recursive formula 
to calculate the total two point impedance was derived. The impedance matrix method was 
applied to obtain the total surface impedance matrix calculated recursively, layer by layer 
by employing the recursive formula N — 1 times for an N layered system. In the case of 
more complex inhomogeneity and a large number of cylindrical layers, the recursive algorithm 
and the alternative integration technique are far superior to methods involving finite element 
analysis and can be performed with completely anisotropic materials. 

Application of the impedance matrix simplifies the formulation of various scattering and 
boundary value problems for cylindrical structures. The impedance matrix can be applied to 
solve various acoustic and elastic scattering problems for arbitrarily layered cylindrical shells 
and solids. As an example acoustic scattering from a solid aluminium cylinder immersed in 
water was considered using the impedance matrix and compared with the literature. Plots of 
the total pressure field, form function and total scattering cross section agree with previously 
published results. 



The equilibrium equations of linear elastodynamics in cylindrical coordinates are [16] 



where u = (u r , ug, u z ) T is the displacement vector, p is the mass density, t r , t^, t z , are 
the traction vectors, K is the 3x3 matrix with K12 = — 1, K21 = 1 and zero otherwise, T 
denotes transpose and a comma suffix denotes partial differentiation. Using the index notation 
(r, 9, z) <f-> (1, 2, 3) the stress strain relationship is 



Appendix A. Cylindrically anisotropic media 



(rt r ) tr + t e ,e + Kt + rt Z)Z = rpii, 



(42) 




(43) 



where the elastic stiffness tensor elements have the usual symmetries C^u = Cjiki = Ckuj, 
and the notation of summation over repeated indices is assumed. Combining the displacement 
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strain equations in cylindrical coordinates with Eq. (I43p the traction vectors are [16] 





(44) 



The form of the various matrices in (JUJ) are, using Voigt notation for the moduli, 



Q = Cie C m C 56 , f = C 26 C 22 C 2 4 I , M = I C 45 C 44 C u , (45) 



C u 


Cl6 


Cl5 


Cl6 




C*56 


C*15 


C56 


C55 


Cl6 


Cl2 


C14 


^66 


C26 


C*46 


C56 


C*25 


C*45 



C*66 


C*26 


C*46 


C*26 


C*22 


C*24 


C46 


C24 


C44 


Cl5 


C14 


Cl3 




C*46 


C36 


C*55 


C*45 


C*35 



r C55 


C45 


C35 


C45 


C44 


C34 


VC35 


C34 


C33 


C56 


C46 


C36 ' 


C25 


C24 


C*23 


C45 


C44 





R = 

We consider cylindrically anisotropic media, illustrated in Figure (TJ for which the density, 
p, and elasticity tensor, C, depend on the radial coordinate r. Solutions are sought in the 
form of time-harmonic cylindrical waves where the displacement and radial traction vectors 
are of the form 

u = U (n) (r ) e i (^+ fc ^-^) , irt r = V (n) ^y^e+k^t) ^ ^ 

where u is the frequency, k z is the axial wave number, and n — 0, 1, 2, ... is the circumferential 
number. Plugging (146j) into (1421) and (14*41 yields the state space representation [16] for the 
state vector rj^ consisting of the displacement and radial traction vectors, 

A^W( r ) = I G (r)^(r) with »? (B) (r) = (v(»)(rj) ■ ^ 
The superscript (n) is omitted for the remainder of the paper. The 6x6 system matrix is 

m 

iG =(i S g{ 3} _g{i}+)> m 

where all terms depend on the radial coordinate, r, and superscript + denotes the adjoint or 
conjugate transpose of a matrix. The 3x3 matrices in f|4*8|) are 

g« = -Q _1 R - i^rQ^P, g< 2 > = = -Q-\ 

g {3} = g {3}+ = f_ r +q i R + ikzT [p^Q^iR _ s - (P r Q *R - S) + ] (49) 
+ r 2 [^(M-P r Q" 1 P)-p W 2 l], 
I is the 3x3 identity matrix, + indicates Hermitian transpose, and 

R = Rk, S = kS, T = T + = k + Tk, with k = K + ml ( = (50) 
The system matrix, G has the important symmetry [16] which follows from the form of (jUJ) , 
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The problem is now reduced to finding a solution to Eq. (j4"7|) i subject to appropriate boundary 
conditions. Next we introduce the matricant and impedance matrices based on solutions of 



Appendix B. Impedance for uniform transverse isotropy 

We consider transversely isotropic solids with the symmetry axis in the z— direction. The 
displacement vector may be decomposed using Buchwald's scalar potentials [19], the functions 

<P, X and ip, 

u = V V + V x ( X e z ) + e z . (52) 

The general solution of the equilibrium equations for transverse isotropy are of the form 

{ip, x, ip} = {fi, X-, ?/;}e l(ri£,+fczZ ~ wt) where 

V = RL ^fn(hr) + ^S l on fi(k 2 r), x 

K\ K 2 - rpl 



X 



TLjrfn(hr), (53) 



^=^R l on fn(hr) + S l on ^ fakir), " ^ 

K\ k 2 

and R l on) S l on) T G ' n are unknown coefficients, fn(x) are cylindrical functions: f^{x) = J n {x) 
for solutions that are regular at r = 0, f 2 (x) = Y n (x) for real valued irregular solutions 
at r = 0, fn{x) = Hn\x) for outgoing (radiating) solutions, f£{ x ) — Hn\x) for ingoing 
solutions, where J n (x) - Bessel function of the first kind; Y n (x) - Bessel function of the second 
kind; Hn' 2 \x) - Hankel functions of the first and second kind. The displacement field can 
be represented as a linear combination of any two of the four types of cylindrical functions 
fn(x), (I = 1,4). The wavenumbers k%, k 2 , k 3 , and non-dimensional numbers k±, k 2 are given 
by 



^1,2 — 7, > ^3 ~~ ~i K i ~ ~7 ~. VT^"' (* ~~ 1 ' 2) , 

2 Cn C44 C 6 6 (C13 + C44 ) k z 

a = (cu + c 44 )pw 2 + (c^ 3 + 2ci 3 c 44 - cn c 33 )k 2 , 
b = 4cn c 44 (puj 2 - c 33 k 2 z ) (pu 2 - c 44 k 2 ) . 

For isotropic material wavenumbers ki, Ki reduce to k\ = u 2 p/(X + 2/i) — k 2 , k\ = k 2 = 
u} 2 p/n — k 2 z1 K\ — 1, k 2 = —k 2 /k 2 . 

The displacement and traction vectors U and V of (j46p are obtained in matrix form for 
each n as 

(Ron\ 
SL , (55) 



a 



=F Va 2 — b 2 _ P^ 2 — c^\k 2 _ C66 k\ — c\\ k. 
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where the summation on / is over any two of the possible 1 = 1,4, and 



X l (r) 



ti'(hr) ft' far) 



k\r J n 
k 



ffL(hr) 
iz l (r)X l (r) 



fHihr) 




(56) 
(57) 



and z , / = 1,4, follows from [17] : 

z l {r) = 



+ c 



2c66 \n2c§Q i/c 2 rc44 
m2c m 2c 66 
i/c 2 rc44 Z z 

6(?/i - 2/2) m(yi - y 2 ) 
-m(yi - y 2 ) bm ~ 62/2 
-16(6 - 6) "(6 - 6) 



i6(£i - 6)' 
«(6 - 6) 




(58) 



C44 



n 2 {iiyi - 62/2) - 61665(3/1 - 1/2) 



6(62/1-612/2) -n 2 (yi 



Co 



c 66 /c 2 r 2 



2/2) 



Vi 

3 m f ) 



1,2), 



(59) 



(J = 1,2,3). 



6(62/1 - 62/2) - ™ 2 (2/i - 2/2) : 
The formula for X z follows by substituting the potentials ( 15B"|) into Eq. (152"|) . The derivation 
of the matrix z l (r) can be found in [17]. Note that z x (r) (/ = 1) is the exact form of the 
conditional impedance of a solid cylinder, i.e. with material at r = and hence bounded 
displacements there [T7] . 

The explicit form of the two point impedance matrix (see Eq. ([H])) of a given transversely 
isotropic layer is 



rk,r k -ij 



z i 

Z k 3 



Z\ 

Z 4 



-Y 1 (r fc _i^ 



-Y 3 (rfc_] 
Y 3 (r fc ) 



X^rfc-i) X 3 (r fc _x 
XHr fc ) X 3 (r fc ) 



(60) 



Equation ( 1601) . which defines the impedance matrix Z, is similar to Eq. (7) of [9] (for the 
stiffness matrix K), and the first and the second matrices on the right hand side of Eq. fl60l) 
are similar to the matrices and (E^J -1 in [91 Eqs. (5) and (3)]. One reason why we prefer 
to use the impedance matrix Z rather than the stiffness matrix as in [9] is that the impedance 
is always Hermitian: Z = Z + . 
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